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Abstract 

We study the spin-1/2 Ising model and the Blume-Capel model at various values of the 
parameter D on the simple cubic lattice. To this end we perform Monte Carlo simulations using 
a hybrid of the local Metropolis, the single cluster and the wall cluster algorithm. Using finite 
size scaling we determine the value D* = 0.656(20) of the parameter D, where leading corrections 
to scaling vanish. We find uj = 0.832(6) for the exponent of leading corrections to scaling. In 
order to compute accurate estimates of critical exponents we construct improved observables 
that have a small amplitude of the leading correction for any model. Analyzing data obtained 
for D = 0.641 and 0.655 on lattices of a linear size up to L = 360 we obtain v = 0.63002(10) 
and r] = 0.03627(10). We compare our results with those obtained from previous Monte Carlo 
simulations and high temperature series expansions of lattice models, by using field theoretic 
methods and experiments. 
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I. INTRODUCTION 



In the neighborhood of a second order phase transition various quantities diverge. For 
example the correlation length, which characterizes the decay of the two-point correlation 
function, behaves as 



e = f±\t\"' X (^1 + b±\tf + ct + d±\tf + e±\tr + ...) , (1) 

where t = (T — Tc)/Tc is the reduced temperature, /+ and /_ are the amplitudes in the 
high and the low temperature phase, respectively, and u is the critical exponent of the 
correlation length. These power laws are affected by confluent corrections, such as 6±|t|^, 
(^±|^r 5 c±|^P^; and non-confiuent ones such as ct. Critical exponents such as u and ratios 
of amplitudes such as f+/f- are universal. This means that they assume exactly the same 
value for any system within a given universality class. Also correction exponents such as 
6 = uu ~ 0.5 and ratios of correction amplitudes as b+/b- are universal. A universality 
class is characterized by the dimension of the system, the range of the interaction and 
the symmetry of the order parameter. The critical exponents a of the specific heat, 7 of 
the magnetic susceptibihty, 77 of the two-point correlation function at the critical point, u 
of the correlation length, S of the magnetization at the critical temperature as a function 
of the external field and (3 of the spontaneous magnetization at a vanishing external field 
are related by so called scaling and hyperscaling relations. This allows to deduce all 
of them from two independent exponents. For reviews on critical phenomena and the 
Renormalization Group (RG) see e.g. 

Here we are concerned with three dimensions, short range interactions and a Z2 sym- 
metry of the order parameter. The best know model undergoing a phase transition in 
this universality class is the spin-1/2 Ising model in three dimensions with nearest neigh- 
bor interactions. Therefore this universality class is called the three-dimensional Ising 
universality class. This universality class is supposed to be realized in a huge range of 
experimental systems: binary mixtures, uniaxial magnets or micellar systems; see e.g. 
|^-l6|. Typically the estimates of critical exponents extracted from experimental data 
are less accurate than those obtained by using the theoretical methods discussed below. 
For example, recent experimental estimates obtained from turbidity data for a methanol- 
cyclohexane mixture are u = 0.632(2) and t] = 0.041(5) Q. 

Critical exponents and amplitude ratios have been computed by various theoretical 
methods such as field theoretic methods or high temperature series expansions and Monte 
Carlo simulations of lattice models. First let us briefly discuss results obtained by the 
e-expansion [8], where the dimension d of the system is given hj d = 4 — e and the 
perturbative expansion in d = 3 joj. The e-expansion of critical exponents has been 
computed up to 0{ef ) [lo[ while the perturbative expansion in d = 3 has been computed 
up to seven loops for the Ising universality class. Since both expansions are divergent, 
some kind of resummation is needed to extract numerical results for critical exponents. In 
the case of the e-expansion the estimates reported in the literature are consistent among 



each other. As a representative result we report in table [T] the one of ref . 12| . In table 
[I we also give results obtained from the perturbative expansion in d = 3 using different 
resummation techniques. For a more complete compilation see e.g. ref. In table |T] we 
give the exponents u, rj and the correction exponent u, since these are directly computed 
by using field theoretic methods. In addition we report the value of 7 that can be compared 
with the results of the high temperature series expansions reported below. Typically, the 
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TABLE I. Numerical results for the critical exponents v, 7, r] and uj obtained by using field 
theoretic methods. The list is by far not exhaustive. We try to give extreme examples; both 
concerning the values found as well as the quoted error bar. In the case of the e-expansion we 
have taken the results that fulfil the boundary condition that for e = 2 the correct 2D Ising 
results are obtained. 



ref year Method u 'j r] oj 



[12] 1998 e-exp 0.6305(25) 

yj 1991 3D exp 0.630 

[11] 1991 3D exp 0.6301(5) 

[12] 1998 3D exp 0.6304(13) 

[14] 1999 3D exp 0.6305 

[15] 2001 3D exp 0.6303(8) 

[16] 2008 3D exp 0.6306(5) 



1.2380(50) 0.0365(50) 0.814(18) 

1.238 0.0355 0.845 
1.2378(6) 0.0355(9) 

1.2396(13) 0.0335(25) 0.799(11) 

1.241 0.0347(1) 0.805 

1.2403(8) 0.0335(6) 0.792(3) 

1.2411(6) 0.0318(3) 0.782(5) 



errors reported for the critical exponents obtained from the perturbative expansion in 
d = 3 are smaller than those obtained from the e-expansion. While the estimates for u 
are all consistent within the quoted errors, clear variations can be observed for t], 7 and 
u. For a discussion of the different resummation schemes that have been used, we refer 



the reader to ref. 16 



In table [TTl we summarize recent results obtained from lattice models. For an exhaustive 
summary of previous works see ref. [3]. The authors of [13] have analyzed the high 
temperature series expansion of improved models on the simple cubic lattice up to 0(/3^^), 
where /3 = l/ksT is the inverse temperature. One of these models is studied here using 
Monte Carlo simulations. Improved means that the amplitudes of leading corrections to 
scaling such as b± in eq. ([T]) vanish. The authors of [18] have studied the high temperature 
series expansion of spin-S Ising models on the simple cubic and the body centered cubic 
lattice up to 0(/3^^). Note that in the spin-S Ising model the spin-variable might assume 
the values —S, —S + 1,...,S' — 1, S. In ref. [l9| the same authors have studied the 

model on the simple cubic and the body centered cubic lattice also up to 0(/3^^). 
These results from high temperature series expansions are all compatible among each 
other. Note that these expansions were performed for lattice models with quite different 
Hamiltonians. Furthermore, there are results for both simple cubic and body centered 
cubic lattices. It is highly plausible that corrections to scaling have different amplitudes 
in these different models. Therefore the agreement of the results gives us confidence 
that there are no undetected systematic errors due to leading, or in the case of improved 
models, subleading corrections to scaling. The results for the exponents u, 7 and t] 
obtained from the high temperature series expansion are clearly more precise than those 
obtained by using field theoretic methods. The results obtained for u using field theoretic 
methods and high temperature series expansions of lattice models are consistent. In the 
case of 7 and rj some of the results obtained by resumming the perturbative expansion 
in three dimensions can be clearly ruled out by the high temperature series expansion. 
Unfortunately, the analysis of the high temperature series expansions does not provide an 
accurate estimate for the correction exponent u. 

Lattice models can also be studied by using Monte Carlo simulations. The finite size 
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TABLE II. Numerical results for the critical exponents u, 7, r/ and uj obtained by analyzing 
high temperature series (HT) and Monte Carlo (MC) simulations of lattice models in the Ising 
universality class. In the case of the Monte Carlo (MC) simulations, some of the authors have 
quoted the statistical and the systematical errors of v and r\ separately. The numbers marked 
by * are not directly given by the authors but are computed by using the scaling relation 
7 = 1^(2 — r]). Note that the error of 7 is computed naively, assuming that the errors of v and 
rj are purely statistical and that the estimates of v and r\ are uncorrelated. For an exhaustive 
summary of previous work see ref. Q/. 



ref 


year 


Method 


I' 


7 


V 


UJ 


[IT] 2002 


HT 


0.63012(16) 


1.2373(2) 


0.03639(15) 


0.825(50) 


[18] 


2002 


HT 


0.6299(2) 


1.2371(1) 


0.0360(8)* 




[19] 


2005 


HT 


0.6301(2) 


1.2373(2) 


0.0363(9)* 




[23] 


1999 


MC 


0.6294(5) [5] 


1.2353(21)* 


0.0374(6) [6] 


0.87(9) 


m 


1999 


MC 


0.6298(2) [3] 


1.2365(11)* 


0.0366(6) [2] 




[25] 


1999 


MC 


0.6296(3) [4] 


1.2367(15)* 


0.0358(4) [5] 


0.845(10) 


[26] 


1999 


MC 


0.63032(56) 


1.2372(13)* 


0.0372(10) 


0.82(3) 


]27] 


2003 


MC 


0.63020(12) 


1.2372(4)* 


0.0368(2) 


0.821(5) 



scaling (FSS) approach [20-23] is well suited to locate the critical temperature and to 
compute critical exponents. Typically one simulates the model directly at the critical 
point. The critical exponents are then extracted from the scaling of the observables with 
the lattice size. For example, at the critical temperature the magnetic susceptibility 
behaves as 

X = aL^-"^ X (1 + bL-'^ + cL''^' + rfL'^'^ + ...) +B , (2) 

where B is an analytic background and L the linear size of a cubic lattice with periodic 
boundary conditions. An exhaustive summary of previous works is given in table 5 of ref. 

In table [TTl we only quote recent works. In 1999 four finite size scaling studies of lattices 
models in the Ising universality class had been published. The results of these works are 
consistent among each other and the accuracy that had been achieved is similar to that of 
the field theoretic calculations. I like to mention that in [26] a special purpose computer 



for the cluster simulation of the Ising model had been used. In the most recent work [27[ , 
which provides the most accurate estimates so far, 11 different models were studied on 
lattices up to a linear site of L = 128. The results obtained for u, 7 and t] are essentially 
consistent with those obtained from the high temperature series expansions. The estimate 
obtained for u is more accurate than that of the high temperature series expansion and 
it is clearly larger than most of the estimates obtained from the perturbative expansion 
in three dimensions. 

The purpose of the present work is to corroborate the lattice results discussed above 
To this end we shall simulate lattices that are considerably larger than those of ref. 127 



Furthermore we shall use improved observables that have been applied in ref. |28| to 
study Ising models with quenched dilution. Here, improved means that the amplitude of 
the leading correction vanishes for any model. Since these observables are constructed 
numerically, in practice some residual amplitude remains. Using these improved observ- 
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ables in the study of improved models, leading corrections are highly suppressed, allowing 
us to ignore them in the finite size scaling analysis. 

Accurate numerical estimates of critical exponents might serve as benchmark for future 
experiments, see e.g. j6|, the analysis of the perturbative expansion in three dimensions, 
as discussed above, or new theoretical approaches such as new ideas in the so called exact 
renormalization group 29| or the Kallen-Lehmann approach [30i] . 

The outline of the paper is the following. First we define the model and the observ- 
ables that are studied. Then we discuss the Monte Carlo algorithm that has been used. 
We give the details of our numerical study. We estimate the fixed point values of the 
phenomenological couplings and the inverse transition temperatures. We give a numer- 
ical estimate of the correction exponent u and obtain a new estimate of D*, the value 
of the parameter where leading corrections to scaling vanish. Next we construct various 
improved observables. Based on this we compute estimates for the critical exponents u 
and rj. 



II. THE MODEL 

The spin- 1/2 Ising model is characterized by the reduced Hamiltonian 

H = -/3 ^ s^Sy - h^s^ , (3) 

<xy> X 

where the spin might assume the values Sx G {—1, 1}. x = {xo,xi,X2) denotes a site of 
the simple cubic lattice, where Xi G {0, 1, 2, Li — 1}. < xy > denotes a pair of nearest 
neighbors on the lattice. We employ periodic boundary conditions in all directions of the 
lattice. Throughout we shall consider Lq = Li = L2 = L and a vanishing external field 
h = 0. The partition function is given by 

Z = J2^M~H) , (4) 

where Yl{s^} denotes the sum over all configurations. 

The Blume-Capel model is characterized by the reduced Hamiltonian 

= -/3 ^ s,Sy + Z)^4-/i^s, , (5) 

<xy> X X 

where now the spin might assume the values G { — 1,0,1}. In the limit D — )■ — 00 
the "state" s = is completely suppressed, compared with s = ±1, and therefore the 
spin-1/2 Ising model is recovered. In > 2 dimensions the model undergoes a continuous 
phase transition for —oo<D< Dtn at a /3c that depends on D. For D > Dtri the model 



undergoes a first order phase transition. Refs. [31n33| give for the three-dimensional 
simple cubic lattice Dtri ~ 2.006, Dtn ~ 2.05 and Dtri = 2.0313(4), respectively. 

Numerically it has been shown that on the line of second order phase transitions there 
is a point, where leading corrections to scaling vanish. In the following we shall call the 
model at this point "improved model". In ref. [34| we find D* = 0.641(8). One should note 
that no effort was made to estimate the systematical error due to subleading corrections 



to scaling. The authors of [27|, |35| have simulated the model at D = \n2 = 0.693147... . 
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At this value of D corrections to scaling are still small compared with the spin-1/2 Ising 
model. At D = ln2 the Blume-Capel model can be mapped into a spin-1/2 Ising model 
with twice the number of sites. This model can be simulated with a cluster algorithm 
without additional local updates as it is the case for general values of D. 



III. THE OBSERVABLES 



The energy of a given spin configuration is defined as 

E = 



(6) 



<xy> 



This definition is convenient for our purpose. One should note however that it deviates 
from the standard textbook definition. The magnetic susceptibility x ^^^^ the second 
moment correlation length C,2nd are defined as 



X 




(7) 



where V = and 



where 



4sin^ Ti/L 



St. 



(9) 



is the Fourier transform of the correlation function at the lowest non-zero momentum. In 
our simulations, we have measured F for the three directions k = 0, 1, 2 and have averaged 
these three results. 

In addition to elementary quantities like the energy, the magnetization, the specific 
heat or the magnetic susceptibility, we compute a number of so-called phenomenological 
couplings, that means quantities that, in the critical limit, are invariant under RG trans- 
formations. We consider the Binder parameter f/4 and its sixth-order generalization Uq, 
defined as 

- 1^, (10) 

where m = ^ Sx is the magnetization of a given spin configuration. We also consider 
the ratio Rz = Za/Zp of the partition function Za of a system with anti-periodic boundary 
conditions in one of the three directions and the partition function Zp of a system with 
periodic boundary conditions in all directions. Anti-periodic boundary conditions in the 
zero direction are obtained by replacing SxSy by —SxSy in the Hamiltonian for links {xy) 
that connect the boundaries, i.e., for x = [L — l,a;i,a;2) and y = (0,xj,X2). The ratio 



Za/Zp can be efficiently evaluated using the boundary flip algorithm 



Here we use 



a modified version of the boundary fiip algorithm as discussed in appendix A 2 of ref 
[s^]. In the following we shall refer to the RG- invariant quantities U2j, 
= ^2nd/ L using the symbol R. 



Rz = Za/Zp and 
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In our analysis we need the observables as a function of (3 in some neighborhood of the 
simulation point. To this end we have computed the coefficients of the Taylor expansion of 
the observables up to the third order. For example the first derivative of the expectation 
value {A) of an observable A is given by 

^-^ = {AE)-{A){E) . (11) 



IV. THE SIMULATION ALGORITHM 

Analogous to [38|, we have simulated the Blume-Capel model using a hybrid of local 
updates and cluster updates. The cluster algorithm only changes the sign of spins. There- 
fore, in order to get an ergodic algorithm for the Blume-Capel model with finite local 
Metropolis updates are used that also can change the modulus \sx \ of the spins. Following 
(39I even in the case of the spin-1/2 Ising model such a hybrid of local and cluster updates 
is superior to the cluster algorithm alone. The authors of jsof also found that such a hy- 
brid algorithm is much less susceptible to systematic errors caused by the imperfection of 
pseudo-random numbers than a pure cluster algorithm. Here we have used a hybrid of 
local Metropolis updates that are implemented by usin g th e multispin coding technique 
(4o| . single cluster updates [ilj and wall cluster updates |42|. In the single cluster update, 
the cluster that includes a randomly chosen site is fiipped. In contrast, in the wall cluster 
update all clusters that include sites that are part of a given plane (the "wall") of the 
lattice are fiipped. 

Motivated by the multispin coding implementation of the local update we have simu- 
lated N}yit = 64 copies of the system in parallel. In the first stage of our study, we have 
used a single random number sequence for the local Metropolis updates of these Nhu sys- 
tems. This leads to some degradation of the performance. To diminish this problem, we 
have used a modified sequence of the pseudo random numbers in the second stage of our 
study. Details are given below. In the case of the cluster updates we could not make use 
of the multispin coding technique. Therefore we have updated the systems one by one, 
using different random number sequences for each of the systems. 

Let us discuss the implementation of the local Metropolis algorithm in more detail. We 
have implemented the spin G {—1, 0, 1} using two bits. To this end we write = cTxT^, 
where ax € { — 1, 1} and G {0, 1}. In terms of these new variables the partition function 
becomes 

Z = C^^exp if3 ^ axCTyTxTy - D^Tx\ , (12) 

where D = D — \n2. Note that subtracting In 2 corrects for the double counting of the 
Sx = state. 

In our local updating scheme we performed consecutive updates of and r^. In the 
first step, the proposal is given by a'^ = —cTx- It is accepted with the standard Metropolis 
acceptance probability 



Pace = min 



l,exp -2P(JxTx ^ (TyT, 



y'y 

y.nn.x 



(13) 



7 



where y.nn.x means that y is a nearest neighbor of x. In the second step, the proposal is 
given by = 1 — Tj.. A natural choice for the acceptance is 



mm 



l,exp (2r^ - 1] 



(14) 



y.nn.x 



Instead, for technical reasons we have implemented a two stage acceptance step. For 
D = 0.641 and D = 0.655, where \D\ is small, we have chosen 



Pr 



accA 



and 



mm 



acc,2 



l,exp (3{1 - 2r^)cr^ 



y.nn.x 



CyTy 



mm 



l,exp (2r,-l)D 



(15) 



(16) 



Detailed balance can be easily proven by going through the four cases which are given 
by positive or negative arguments of the exponential function in Eqs. ( !T5l[T6|) . We take 
two uncorrelated random numbers ri and r2 from a uniform distribution in [0, 1]. If both 
Pacc,i > ^1 and Pacc,2 > 1^2 the proposal is accepted. 

In the first stage of our study, we have used the same random number sequence for all 
Nhit = 64 systems that we have simulated in parallel. In the second stage of the study, 
we have used a modified sequence of the random numbers for the acceptance step (IT^ . 
We have used a 64 bit integer random number in addition to the random number r that 
is uniformly distributed in [0,1]. If the i^^ bit of this integer random number is 1 we 
take r itself for the acceptance step of the z*'* system. Otherwise, if the bit is 0, we take 
1 — r instead. This modification considerably reduces the correlation among the Nut = 64 
systems that are simulated in parallel. 

We have compared the performance of this local update with that of a local heat 
bath using a standard implementation. To this end, we have simulated a 16^ lattice at 
(3 = 0.3877218, which is close to (3c as we shall see below. The integrated autocorrelation 
times in units of sweeps of the energy density and the magnetic susceptibility are by 
a factor of about 1.3 larger for the Metropolis update discussed here than for the heat 
bath update. One sweep over Nut = 64 systems in parallel using the multispin coding 
technique takes about 4 times as much CPU-time as one sweep over a single system 
using the standard implementation of the heat-bath update. In order to compare the 
efficiency of the two local updates, we have computed the statistical error of the energy 
and the magnetic susceptibility, taking into account the possible correlation among the 
Nut = 64 systems in the case of the multispin coding implementation. To this end we 
have performed a Jackknife analysis, where we have first averaged the measurements of the 
Nut = 64 systems at a given iteration of the Monte Carlo simulation. Taking the inverse 
of the statistical error squared times the CPU time needed as measure of the efficiency, 
we find a performance gain of a factor of about 10 of the multispin coding implementation 
of the local Metropolis update compared with the standard implementation of the heat 
bath update. 

During the simulation, local Metropolis sweeps, single cluster and wall cluster updates 
are performed in a certain sequence. In the following we denote an elementary building 
block of the sequence by cycle. In the case of our most recent simulations {D = 0.655) 
such a cycle is composed of 
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• 4 X (2 Metropolis sweeps followed by L/16 single cluster updates) 

• 3 Metropolis sweeps 

• one wall-cluster update 

In the case of the wall-cluster update we chose the wall to be perpendicular to the 0, 1 
and 2- axis in three subsequent cycles. The position of the wall along the axis is chosen 
randomly each time. The parameters of the cycle are chosen such that roughly the same 
amount of CPU-time is spent in each of its three components. 

In order to study the performance of the algorithm, we have performed preliminary 
simulations for D = 0.655, where we have determined the autocorrelation function p{t) of 
the magnetic susceptibility and the energy density. The statistics of these runs is 300000 
update-cycles for the lattice sizes L = 16, 32, 64, and 128 and 82000 update-cyclcs for 
L = 256 at /3 = 0.3877218, which is close to our final estimate of (3c- The integrated 
autocorrelation time is given by 

t=i 

where we have chosen t^ax — 6t, selfconsistently. Fitting our results for integrated 
autocorrelation times in units of update-cycles we get 

= 0.70(4) X L°-=^^(i) (18) 

for the magnetic susceptibility and 

te = 0.47(2) X L°-^2(i) ^19^ 

for the energy density. This means that the autocorrelation times are only a few cycles, 
even for our largest lattices. 

We have estimated the statistical errors of the observables using the Jackknife method. 
As input of this analysis we have taken data that are already averaged over the N^it — 64 
systems that are simulated in parallel and might be correlated by the use of a common 
sequence of random numbers during the Metropolis updates. Therefore the possible cor- 
relation among these Nbu = 64 systems does not affect the correctness of the estimate of 
the statistical errors. 

To figure out how much this correlation does affect the efficiency of the algorithm, 

we have computed the statistical error, taking only one system, and for comparison, 
averaging over all Nm systems. If the simulations were independent, the square of the 
ratio of these errors, denoted by in the following, would be equal to iV^jf. In fact we see 
some performance loss due to the use of a common random number sequence. For L = 16 
wc get for the energy density ~ 28.6 and for the magnetic susceptibility R^ ^ 36.2. 
Fortunately these numbers increase with increasing lattice size. For L = 256 we get for 
the energy density R^ ~ 42.2 and R^ ~ 48.8 for the magnetic susceptibility. 

In order to give an accurate result for the performance gain that is achieved by using 
our particular multispin coding implementation of the local update, one would have to 
tune the parameters of the update cycle for both types of the local update. For lack of 
time this could not be done. Since the local update is only one of the three components 
of the complete update cycle, likely the gain is moderate, certainly less than a factor of 
two. 
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A. The simulations: CPU time and statistics 



In a first stage of tlie study we liave simulated tlie spin-1/2 Ising model and tlie Blume- 
Capel model at D = 0.641, In 2, 1.15 and 1.5. Note that D = 0.641 is the estimate of ref. 



for D* and D = ln2 has been simulated before by the authors of refs. [27|, |35j. At 
D = 1.15 the amplitude of leading corrections to scaling has about the same magnitude 
as for the spin-1/2 Ising model but opposite sign. A preliminary analysis of these data 
resulted in D* ?a 0.655. Therefore we have simulated at D = 0.655 in a second stage of 
our study. 

We have simulated lattices of a linear size L up to L^ax = 96, 200, 360, 300, 64 and 48 
for the spin-1/2 Ising model and the Blume-Capel model at D = 0.641, 0.655, In 2, 1.15 
and 1.5, respectively. In table UTTl we have summarized in detail the lattice sizes that we 
have simulated and the statistics of these simulations. 

In total we have spent the equivalent of 3.5, 9, 16, 3, 3, 0.1 CPU years on a single core of 
a Quad-Core AMD Opteron(tm) Processor 2378 running at 2.4 GHz for the spin-1/2 Ising 
model and the Blume-Capel model at D = 0.641, 0.655, In 2, 1.15 and 1.5, respectively. 

As random number generator we have used the SIMD-oriented Fast Mersenne Twister 



algorithm [43|. As a check we have repeated our simulations at -D = 0.655 using the 
WELL Random number generator [3] with about one third of the statistics reported 
in table IIIII In particular we have used the program 'WELL44497a.c" provided by the 
authors. We found that the estimates of individual observables are consistent. We have 
also repeated part of the finite size scaling analysis using these data. For given ansatze we 
found consistent, even though less precise results for the critical exponents. One should 
note that the statistical error of the fit parameters is often much smaller than the final 
error that also includes systematical errors due to subleading corrections. The following 
analysis is only based on the simulations using the Mersenne Twister algorithm [4^ . 



V. f3c AND THE FIXED POINT VALUES OF PHENOMENOLOGICAL COU- 
PLINGS 

In a first step of the analysis we have studied the finite size scaling behavior of the 
phenomenological couplings at D = 0.655, since here we have accumulated the best 
statistics and secondly, as we shall see below, this value of D is closest to D* among the 
values that we have simulated. 

At the critical point a phenomenological coupling behaves as 

R{L, =R* + aL-'^ + bL-^' + cL-^'' + ... , (20) 

where u ~ 0.8 as discussed in the introduction. Below we shall find u = 0.832(6). The 
subleading corrections exponent is u' = 1.67(11) j45|. Furthermore, there should be 
corrections with oo" ~ 2 due to the breaking of the rotational symmetry by the lattice 
or due to the analytic background of the magnetic susceptibihty. Motivated by eq. (l20l) . 
we have fitted our data with three different ansatze 

R{L,f3,)=R* (21) 
R{L, (3c) = R* + aL-'' (22) 
R{L, =R* + aL-'' + bL~'^ , (23) 



10 



TABLE III. We give the number of update-cycles divided by 64 x 15000 as a function of the 
lattice size and the value of the parameter D. For a discussion see the text. 



L 


Ising 


0.641 0.655 


In 2 


1.15 


1.5 


10 


10000 


10000 


4005 




10000 




11 






4005 








12 


9593 


20000 


4000 


4083 


10000 


1000 


13 






4005 








14 


11747 


10000 


4005 


3003 


10000 




15 






3994 








16 


9740 


20200 


3999 


2371 


9917 


1000 


17 






4000 








18 


11524 


10000 


3993 


1807 


11102 
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where we have used in eq. f l22l) the choices ei = 0.83, ei = 1.6 or ei = 2 and in eq. fl23l) 
ei = 0.83 and e2 = 1.6 or €2 = 2. Here and in the following ansatze, we denote a correction 
exponent with a fixed value by e. Instead, if it is a free parameter we shall denote it, as 
usual, by u. Here we need the phenomenological couplings i? as a function of the inverse 
temperature. To this end we have used the Taylor-expansion around the value Ps of the 
inverse temperature that we have used in the simulation. We have checked that the result 
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TABLE IV. Fitting the data for Za/Zp obtained at D = 0.655 with the ansatze m\22\23]) . Lr, 
is the minimal lattice size that is included into the tit. For a discussion see the text. 



ansatz ei £2 Lmin {Zg/Zp)* x /d-o.f- 

ED - - 32 0.387721745(10) 0.542489(14) 9.5/10 

[22] 0.83 - 18 0.387721730(12) 0.542589(33) 12.4/14 

[22] 1.6 - 12 0.387721729(10) 0.542558(12) 24.0/20 

[22] 2 - 10 0.387721734(10) 0.542532(8) 27.4/22 

[23] 0.83 1.6 10 0.387721746(12) 0.542448(46) 27.6/21 

[23] 0.83 2 10 0.387721740(12) 0.542502(38) 26.8/21 



TABLE V. Results for the inverse critical temperature f3c Sit D = 0.655 obtained from the FSS 
study of various phenomenological couplings. In addition we give the fixed point values R* of 
these quantities. For a discussion see the text. 

Zg/Zp ^2nd/ L U% 

Pc 0.38772174(2) 0.38772174(2) 0.38772173(2) 0.38772173(2) 
R* 0.5425(1) 0.6431(1) 1.6036(1) 3.1053(5) 



for /3c and are sufficiently close to avoid significant truncation effects. This way, for 
example eq. ([2T|) becomes 

i?(L,/3,) = /?*-ci(L,/3,)(/3,-/3,)-^^^:|^(/3,-/3,)2-^^^ , (24) 

where R* and /3c are the two parameters of the fit. Since we have chosen as a good 
approximation of /3c, we could ignore the relatively small statistical error of the Taylor 
coefficients ci, C2 and C3, which simplifies the fit. 

As an example let us discuss the results obtained for Zg/Zp in more detail. A selection 
of our results is given in table IIVI We have fitted the data for all linear lattice sizes L 
that are larger than or equal to a certain L^j^. Starting from the Lmin given in column 
5 of table [IV] the x^/d-o.f. is close to one. 

Taking into account the variation of the results over the different ansatze we arrive at 
our final estimate /3c = 0.38772174(2) and {Zg/Zp)* = 0.5425(1). We performed a similar 
analysis for C,2nd/L, U4 and Uq. Our final results are summarized in table [V] We find that 
the estimates of /3c obtained from different phenomenological couplings are consistent 
within error bars. We take the average 

/3c = 0.387721735(25) (25) 

as our final estimate of the inverse critical temperature. The error bar is chosen such that 
it covers all results given in table [V] including their error bars. 

Our result for is about 3 times the combined error smaller than = 1/0.62341(4) = 



1.60408(10) given in [27[. We regard our result as more reliable, since we have simulated 
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TABLE VI. Estimates of the inverse critical temperature (3c of the spin-1/2 Ising model and the 
Blume-Capel model at various values of D. For a discussion see the text. 



D Pc 

Ising 0.22165463(8) 
0.641 0.38567122(5) 
0.655 0.387721735(25) 
In 2 = 0.69314718... 0.39342239(8) 
1.15 0.4756110(2) 
1.5 0.5575303(10) 



larger lattices and have carefully estimated systematic errors due to subleading corrections 
that are not included into the fit. 

In the case of the other models we also determined /3c by fitting with the ansatze f l2T]|22|23p 
Here however we have used the results for (Z^/Zp)*, {^2nd/L)*, Ul and f/g given in ta- 
ble |V] as input. The final results obtained this way are summarized in table IVII For 
completeness we have included the results for D = 0.655 given in eq. fl25l) . 

Our result for /3c of the spin-1/2 Ising model is fully consistent with Pc = 0.22165455(3) 
given in [23|. For a sumr nary of previous results for /3c of the spin-1/2 Ising model we refer 
the reader to table 1 of [34]. Our result for /3 c at D = In 2 is by 1.5 times the combined 
error larger than /3c = 0.39342225(5) given in 27 . 



VI. THE CORRECTION EXPONENT uj AND THE IMPROVED MODEL 

In this section we study the cumulants t/4 and Uq at a fixed value of Za/Zp or ^2nd/L. 
To this end one determines the inverse temperature Pf{L) defined by 

R,i^PJiL)) = R^J , (26) 

where -Ri is either Za/Zp or ^2nd/L and Rj^i the required value. As we take the fixed 
point values of Za/Zp and ^^nd/L obtained above. We define 

R2{L)=R2{L,f3f{L)) , (27) 

where R2 is, in our case, either U4, or f/g. In the following we shall denote R2 by R2 at 
Ri = Rij. In practice we have done these calculations using the Taylor-expansion of Ri 
and R2 around the value Ps that we have used in the simulation up to third order. We 
have checked carefully that Ps and /3j are sufficiently close to avoid significant truncation 
errors. 

One finds, see e.g. section III of js^] 

R{D, L)=R* + a{D)L-'' + b{D)L-'^' + ... +c a\D)L-^'' + ... , (28) 

where we should note that the correction amplitudes depend on the parameter D of 
our model. The improved model is characterized by a vanishing amplitude of leading 
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corrections to scaling. Hence D* is given by the zero of a{D). We have analyzed the data 
of 5 different models in combined fits: The Ising model and the Blume-Capel model at 
D = 0.641, D = 0.655, D = ln2 and D = 1.15. To this end we have employed various 
ansatze that are derived from eq. ( 128|) : 

R{D,L) = R* + a{D)L-'^ (29) 

R{D, L) = R* + a{D)L~'^ + c a^{D)L-'^^ (30) 

R{D, L) = R* + a{D)L-'^ + c a^{D)L-^'^ + hL'^ (31) 

R{D, L) = R* + a{D)L-'^ + c a^{D)L-'^'^ + d a^{D)L-^'^ (32) 

R{D, L) = R* + a{D)L-'^ + c a^{D)L-^'^ + d a^{D)L-'"^ + bL'' . (33) 

In the ansatz fl29|) the free parameters of the fit are R*, a(Ising), a(0.641), a(0.655), a(ln2), 
a(1.15) and the correction exponent u. In the ansatz f l5U]) we have in addition the param- 
eter c. In the ansatz flHT]) we have added the term bL^"^ to take subleading corrections into 
account. Here we make the approximation that the parameter b is model independent. 
We fix the subleading correction exponent e = 1.6 or e = 2. In the ansatz ([32]) we take 
into account corrections oc L~^'^ . Finally in the ansatz ( 133|) we add, similar to eq. ( 13T|) a 
term bL~'^. 

In the case of the ansatz f l29p fits with x^/d.o.f.< 2 are only obtained for Lmin > 36. 
Instead, fitting with ansatz (ED]) we get for at Za/Zp = 0.5425 xVd-o.f.= 62.4/62 
already for Lmin = 16. The results for the parameters of this fit are u = 0.832(1), 

= 1.60357(1), a(Ising) = -0.2983(6), a(0.641) = -0.0067(2), a(0.655) = -0.0006(2), 
a(ln2) = 0.0167(2), a(1.15) = 0.380(1), and c = 2.08(3). Note that here and in the fol- 
lowing the errors quoted for results of individual fits are purely statistical. Extrapolating 
a(0.641) and a(0.655) we get D* = 0.6564(5). 

We estimated the systematic error due to corrections that are not taken into account 
in the ansatz (!30!) from the variation of the results obtained with the ansatze fl31ll32f33l) 
and by using Uq instead of 11^. Furthermore we have redone the analysis for Ui and Uq 
at ^2nd/L = 0.6431. We arrive at the final estimates 

tu = 0.832(6) (34) 
D* = 0.656(20) . (35) 

It also follows from the fits that the amplitude of corrections to scaling at D = 0.655 is 
at least by a factor of 30 smaller than that of the spin-1/2 Ising model. 

VII. IMPROVED OBSERVABLES 

The exponent u can be obtained from the behavior of the slope of a phenomenological 
coupling at the critical point: 



OR 

'dp 



= aL^''' (l + 6L-" + ...) . (36) 



The exponent rj can be extracted from the behavior of the magnetic susceptibility at the 
critical point: 

x\p=p^ = aL^-'^ (1 + 6L-- + ...) . (37) 
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Note that the coefficients a and b of course take different values in eq. (15^ and eq. 
Such a procedure requires an estimate of f3c- To avoid this we have studied, following [25| 
the slopes and the magnetic susceptibility at f3f as defined in eq. f l26|) . These quantities 
behave as 

= a{D)L^''' {l + h{D)L-^ + ...) (38) 



dR 



dp 



and 



X 



(39) 



X\0=i3f 

Again we have computed these quantities using their Taylor-expansion around Ps up to 
the third order. 

Here, following j28|, we shall study improved versions of the slopes and the magnetic 
susceptibility. This means in the ideal case that the amplitude of leading corrections 
vanishes for any model. In practice, as we shall see below, we can construct quantities 
for that the amplitude of leading corrections is suppressed by more than one order of 
magnitude. Using such quantities in the case of improved models ensures that leading 
corrections to scaling are suppressed by two to three orders of magnitude compared with 
standard observables in the case of e.g. the spin-1/2 Ising model. This is sufficient to 
ignore leading corrections to scaling in the analysis of our data. 

Let us discuss in detail the construction of the improved observable at the example of 
the magnetic susceptibility. We consider 



(40) 



where x is chosen such that the amplitude of leading corrections vanishes. Note that 
instead of O4 also Uq could be used. It is important to take a phenomenological coupling, 
where leading corrections to scaling are clearly visible. Let us recall the finite size scaling 
behavior of the Binder cumulant 



Ui{L,D) = U: + buiD)L-'^ + ... . 
Inserting eq. (15^ and eq. fllTl) into eq. fjlOl) we get 

buiD) 



X^mp{L,D) = a{D)U^L 



X T 2—ri 



1 + 



b{D)+x- 



L-^ + ... 



(41) 



(42) 



Hence the exponent defining the improved observable is given by 



X 



TT* 

-b{D) 



buiD) 



(43) 



Note that ratios of correction amplitudes are universal. Therefore the exponent x does 
not depend on D. It can be best determined by analyzing data obtained for models with 
relatively large corrections to scaling. For example one might consider the spin-1/2 Ising 
model to this end. We have already determined 6[/ (Ising) and U^* in the previous section. 
In order to obtain b{D) one would fit x{L, D) with ansatze motivated by eq. (139|) . 

However it turns out to be more efficient to study ratios of observables taken from 
two different models. This way, critical exponents cancel and therefore fits have less 
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where now 



parameters and become more rehable. In particular we shall study the spin-1/2 Ising 
model and the Blume-Capel at D = 1.15. We define 

and 

^ ^ U{L,D = 1.15) [/* ' ^ ^ 

X = -[6(Ising) - 6(1.15)] . ^\ . (46) 

6;7(Ismg) - 6c/(l-15) 

The exponent x can be directly obtained from fits with the ansatz 

Ru{LYR^{L) = C , (47) 

where x and C are the parameters of the fit. To check for the effect of sub leading 
corrections we have also fitted the data with the ansatz 

Ru{LfR^{L) = C + cL~' , (48) 

where c is an additional parameter and we have fixed either e = 1.6 or e = 2. Fixing 
Za/Zp = 0.5425, fits with the ansatz P7|) have an x^/d.o.f. ~ 1 starting with Lmm = 16. 
Using Lmin = 16 we get x = —0.656(1). Instead using the ansatz (HHj) we get x^/d-o.f. ~ 1 
already for L^m = 10. The results for L^m = 10 are x = —0.665(2) and x = —0.661(2) 
for e = 1.6 and e = 2, respectively. As our final result we quote x = —0.66(1), where the 
error is chosen such that it covers the three estimates given above. In a similar fashion 
we arrive at a; = —0.57(2) for fixing ^2nd/ L = 0.6431. 

In figure [1] we demonstrate the effectiveness of the improvement. We have analyzed 
our data for x at Za/Zp = 0.5425 for the Ising model and the Blume-Capel model at 
D = 1.15. To this end, we have fitted our data with the ansatz 

X = aL^-^ + B , (49) 

where B is an analytic background. Using the standard magnetic susceptibility, we get 
X^/d.o.f.= 4.6/4 for Lmin = 32 in the case of the Ising model and x^/d-o.f.= 2.6/5 for 
Lmin = 24 in the case of the Blume-Capel model at D = 1.15. Nevertheless for e.g. 
Lmin = 32 the results for t] obtained from the two different models differ by more than 
20 standard deviations. In contrast, for the improved magnetic susceptibility the results 
obtained for the two models are quite similar. In particular for Lmin = 24 the estimates for 
7] obtained from the Ising model and the Blume-Capel model at D = 1.15 are consistent 
within the error bars. 

We also have constructed improved slopes 

Simp{L, D) = U,{L, DrS{L, D) , (50) 

where x is chosen such that leading corrections to scaling vanish. We have determined x 
analogous to the case of the magnetic susceptibility discussed above. To this end we have 
computed the ratios 

g(L, Ising) 
""^^"^^ = ^(L,D = 1.15) • 
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FIG. 1. Results for the critical exponent rj obtained by fitting the standard and the improved 
magnetic susceptibility at Za/Zp = 0.5425 for the Ising model and the Blume-Capel model at 
D = 1.15 using the ansatz (|49p . L^in is the minimal lattice size that is taken into account. 
In the case of the improved magnetic susceptibility, the results obtained from the two different 
models fall nicely on top of each other. The dashed lines should only guide the eye. 



TABLE VII. Exponent x of improved slopes as defined by eq. (I50|) . In the first column we give 
the phenomenological coupling and its value that is used to dehne Pj. In the first row we give 
the quantity whose slope is considered. For a discussion see the text. 



fix ; slope of Zg/Zp C2nd/L U4 Uq 
Za/Zp = 0.5425 0.52(2) 0.77(3) -1.21(5) -2.73(5) 
^2nd/L = 0.6431 0.54(2) 0.81(2) -1.21(3) -2.71(4) 



As discussed above for the case of the magnetic susceptibility, we have fitted 

Ru{LYRs{L) = C (52) 

with X and C as free parameters and, as check 

Ru{LfRs{L) = C + cL-' , (53) 

where c is an additional parameter and e is fixed to either 1.6 or 2. Our final results for 
the exponent x are summarized in table I VI II 
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TABLE VIII. Exponent x of improved slopes as defined by eq. i54\) . In the first column we give 
the phenomenological coupling and its value that is used to define j3f. In the first row we give 
the quantity whose slope is mixed with that ofU^. For a discussion see the text. 

fix ; slope of Zg/Zp C2nd/L 
Za/Zp = 0.5425 0.29(1) 0.39(1) 
C2nd/L = 0.6431 0.31(1) 0.41(1) 



Furthermore we have constructed quantities of the type 

Sij = \Si\ \Sj\ , (54) 

where x is again chosen such that the amphtude of the leading correction vanishes. Here 
we performed fits with the ansatz 

RlRs7 = C , (55) 
where x and C are the parameters of the fit and 

RsA-'' = C + cL-^ (56) 

with the additional parameter c. Also here we have fixed either e = 1.6 or e = 2. In 
practice we have combined the slope of [/4 with the slope of Za/Zp or ^2nd/L. Our results 
for the exponent x are summarized in table IVIIII Notice that the results obtained for 
Za/Zp = 0.5425 and ^2nd/ L = 0.6431 are similar but not identical. 



VIII. THE EXPONENT r] 

We have fitted our data for the improved magnetic susceptibility at Za/Zp = 0.5425 
and C,2nd/L = 0.6431 using the ansatze 

Xin.p = a{D)L^-^ (57) 
Ximp = a{D)L'-'^ + B{D) (58) 
Ximp = aiD)L'-'^ {l + d{D)L-') . (59) 

In the ansatz (|58ll we have taken into account the analytic background B of the magnetic 
susceptibility. Since 1] is small, the parameter B also takes effectively into account other 
corrections that have a correction exponent u" ~ 2 like for example the breaking of the 
rotational symmetry by the lattice. In the ansatz fl59|) we have set e = 1.6. Using this 
ansatz we try to estimate the possible effect of a correction caused by u' = 1.67(11) jisj 
on our estimate of t]. 

We have fitted the data for D = 0.641 and D = 0.655 in a common fit. The parameters 
of these fits are a(0.641), a(0.655) and rj in the case of the ansatz (|57|) and in addition 
5(0.641) and 5(0.655) or c/(0.641) and c/(0.655) in the case of the ansatz (EH]) or ([59]), 
respectively. In figure |2] we have plotted estimates of t] obtained by fitting with the 
ansatz (1571) as a function of i^mL- Up to Lmin = 48 the results fall roughly on a straight 



18 



0.0365 



0.036 




Za/Zp= 0.5425 
xi2nd/L= 0.6431 



0.0355 



0.035 



I 

0.001 



I 

0.002 



I 

0.003 



0.004 



FIG. 2. Results for the critical exponent ry obtained by fitting the improved magnetic susceptibil- 
ity at Za/Zp = 0.5425 and at i2nd/L = 0.6431 using the ansatz (f57|l . Data for the Blume-Capel 
model at = 0.641 and D = 0.655 are taken into account. Lmin is the minimal lattice size that 
is taken into account. The dashed lines should only guide the eye. For a discussion see the text. 



line, indicating that corrections with an exponent e ~ 2 are present. For Lmin = 128 one 
finds that x^/d-o.f. is smaller than one for both taking the improved susceptibility at 
Za/Zp = 0.5425 and ^2nd/L = 0.6431. The estimate rj = 0.03636(20) covers both results 
obtained at Lmin = 128, including their error bars. 

Next we have fitted our data with the ansatz (!58|) . For both the improved magnetic 
susceptibility at ^2nd = 0.6431 and at Za/Zp = 0.5425 we get x^/d.o.f.< 2 starting from 
Lmin = 14. The results for t] are plotted in figure [31 In the case of ^2nd = 0.6431 the 
estimate of rj is increasing with increasing Lmin, while it is decreasing for Za/Zp = 0.5425. 
For Lmin = 32 the two results are consistent within error bars. 

We read off our final estimate 

77 = 0.03627(10) . (60) 

The error estimate is chosen such that it also covers results obtained with the ansatz fl59l) 
cLnd J^YfitTi — 32. 



IX. THE CRITICAL EXPONENT u 

In order to determine the exponent v we performed combined fits of our data for the 
improved slopes aX D = 0.641 and D = 0.655. In a first step of the analysis we have fitted 
the improved slopes with a power law without any correction 

S = a{D)L^/'' , (61) 
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FIG. 3. Results for the critical exponent ry obtained by fitting the improved magnetic susceptibil- 
ity at Za/Zp = 0.5425 and at i2nd/L = 0.6431 using the ansatz ([58]) . Data for the Blume-Capel 
model at -D = 0.641 and D = 0.655 are taken into account. Lmin is the minimal lattice size that 
is taken into account. The dashed lines should only guide the eye. For a discussion see the text. 



where the amplitudes a(0.641), a(0.655) and the exponent v are the parameters of the fit. 
In figure m we give the results for i/ as a function of i^"^^, where is the minimal lattice 
size that is included into the fit. In the figure we give only results for taking the slopes at 
Za/Zp = 0.5425. Those for the slopes at C,2nd/L = 0.6431 behave in a very similar way. 

We find that the result for u obtained from the improved slope of Za/Zp is increasing 
with increasing Lmin, while the one obtained from ^2nd/L is decreasing. The results 
obtained from the improved slope of [/4 and Uq are quite similar. They only slightly 
increase with increasing Lmin- We have also plotted results obtained from the combined 
slopes (15^ . In the case of combining the slope of f/4 with that of i2nd/L the estimate 
of V is decreasing with increasing Lmin-, while for combining the slope of f/4 with that of 
Za/Zp it is increasing. In all rather large Lmin is need to get acceptable values 

for x^/d.o.f. . In the worst case, for the improved slope of Za/Zp only for Lmin > 56 a 
X^/d.o.f. smaller than two is reached. The behavior of the estimates of v for Lmin < 48 
is consistent with the fact that the dominating corrections have an exponent u' ^ 2. 
For larger Lmin, the variation of our estimates of u with Lmin seems to be dominated by 
statistical fluctuations. For Lmin = 128 we get x^/d.o.f. smaller than one for all quantities 
that we have considered. The estimate u = 0.6301(3) covers the results, including the 



statistical error, of all our fits for Lmin = 128. Note that in ref. [27[ L = 128 is the largest 
lattice size that is simulated. 

Motivated by these observations, we have fitted our data with the ansatz 

S = a{D)L^/'' X (1 + bL~') , (62) 
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FIG. 4. Results for the critical exponent v obtained by fitting improved slopes of various 
phenomenological couplings at Za/Zp = 0.5425 as a function of L^-^, where Lmm is the minimal 
lattice size that is included into the fit. The dashed lines should only guide the eye. For a 
discussion see the text. 



where we have fixed e to either 1.6 or 2. Since already a(0.641) and a(0.655) are very 
similar, we have chosen the parameter b to be model independent. Let us first discuss 
the fits with e = 2. Such fits give x^/d.o.f. close to one already for Lmm = 10. In the 
lower part of figure |5] we have plotted the results obtained from the slopes of different 
quantities for 10 < Lmin < 24. These different estimates of u are consistent among each 
other. Furthermore there is little variation of the results with Lmm- 

In the upper part of figure |5] we plot the corresponding result for e = 1.6. Here the 
X^/d.o.f. is somewhat larger than for e = 2. Also the result for u clearly depends on the 
quantity that is analyzed. We conclude that the numerically dominant corrections have 
an exponent e ~ 2. Motivated by these fits, we take u = 0.63002 as our final result. Since 
we can not exclude that there are also corrections with an exponent e ~ 1.6, we take these 
fits into account in our final error of u. For Lmin = 22 and Lmin = 24 all results that we 
have obtained with the ansatz fl62|) . including their error bar, are contained in the interval 
[0.62992,0.63012]. Therefore we quote as our final result 

= 0.63002(10) . (63) 



X. SUMMARY AND CONCLUSIONS 

We have simulated the spin-1/2 Ising model and the Blume-Capel model on the sim- 
ple cubic lattice using linear lattice sizes L < 360. Using finite size scaling methods we 
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FIG. 5. Results for the critical exponent v obtained by fitting improved slopes of various 
phenomenological couplings at Za/Zp = 0.5425 with the ansatz ()62p as a function of Lmin- In 
the upper part of the figure the correction exponent is fixed to e = 1.6 and in the lower part it 
is fixed to e = 2. The dashed lines should only guide the eye. For a discussion see the text. 
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have determined critical properties of these models. In particular we have determined 
the value D* = 0.656(20) of the parameter D of the Blume-Capel model, where leading 
corrections to scaling vanish. We have accurately determined the inverse of the criti- 
cal temperature for various values of D, in particular /3c(0.641) = 0.38567122(5) and 
^c(0.655) = 0.387721735(25). We have computed the critical exponents u = 0.63002(10) 
and 7] = 0.03627(10) as well as the exponent u = 0.832(6) of leading corrections to scaling. 
The errors quoted for these final results cover statistical as well as systematical errors. 
Systematical errors are due to the fact that power laws like eqs. fl36|37l) that govern the 
finite size scaling behavior of physical quantities at the critical temperature are subject to 
an infinite series of correction terms. Fitting Monte Carlo data, only few of these correc- 
tion terms can be taken into account. In the present study, we have effectively eliminated 
the leading correction oc L~'^ by simulating an improved model and analyzing improved 
observables as discussed in section IVIII In our ansatze we take into account a sub- leading 
correction with the exponent u' = 1.67(11) predicted by j45| or u" ~ 2 due to the break- 



ing of the rotational symmetry by the simple cubic lattice |46j or due to the analytic 
background of the magnetic susceptibility. We estimate the error caused by correction 
terms that are not included by comparing the results obtained by using different ansatze 
and, even more important, by fitting different quantities. One expects that in the generic 
case the amplitudes of corrections are different for different quantities. In the case of the 
critical exponent u we have studied the slope of four different phenomenological couplings: 
The cumulants and Uq, the ratio of partition functions Za/Zp, and the second moment 
correlation length over the linear lattice size C,2nd/L. We regard the estimates of the error 
obtained this way as quite robust and therefore the results obtained here should serve well 
as benchmark for experimental studies as well as new or developing theoretical methods. 

Our results are fully consistent with those obtained from high temperature series ex- 
pansion of lattice models ITMH; See table HI We find a small discrepancy with the 



Monte Carlo results of ref. [27|; See table HIl Note that the authors of l 27| did not take 
into account a sub-leading correction with the exponent u' = 1.67(11) |45| analyzing their 
Monte Carlo data. 

The accuracy that is reached now by lattice methods has clearly outpaced that of field 
theoretic methods. Furthermore, comparing with the numbers that are summarised in 
table [H we notice that most of the results for t] and u obtained from the perturbative 
expansion in three dimensions fixed are at odds with ours^while those of ll|, ll3| are in 
reasonable agreement. Note that, as discussed by Nickel [l3|, the sub leading correction 
exponent u' = 1.67(11) [45] also plays a crucial role in the analysis of the perturbative 
series in three dimensions fixed. Therefore, it would be highly desirable to get an estimate 
of u' by using a different method. 

Using Monte Carlo simulations, the error of the estimates of the critical exponents can 
be further reduced just by spending more CPU time. To this end one has to increase the 
statistics as well as enlarge the size of the lattices that are simulated. Keeping the statis- 
tical error and the systematical one proportional, the effort increases as error"^"'^^"'"^^/'^ 
with a decreasing error, where the first factor error"^ is related to the increased statistics 
and the second to the larger linear lattice size L that is needed to reduce the systematical 
error. Here we assume that the systematical error is proportional to L"'^ , since, as we 
have shown here, leading corrections can be eliminated. The effort at a fixed statistical 
accuracy behaves as L*^"^^, where ci = 3 is the dimension of the system and z is the critical 
dynamical exponent. In a recent study of a spin glass ^] about 1000 years of CPU 
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time on one core of a CPU of similar performance as the one used here had been spent. 
This is about a factor of 30 more CPU time than we have spent here. One should notice 
however that this factor in CPU time only would allow to reduce the errors of the critical 
exponents by a factor of about 2.3, where we have assumed u' ~ 1.6 and z ~ 0.4; see 
Eqs. (1T8][T9|). 
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